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Abstract. We investigate a model of the GnRH pulse and surge generator, with the defi- 
nite aim of constraining the model GnRH output with respect to a physiologically relevant 
list of specifications. The alternating pulse and surge pattern of secretion results from the 
interaction between a GnRH secreting system and a regulating system exhibiting fast-slow 
dynamics. The mechanisms underlying the behavior of the model are reminded from the study 
of the Boundary-Layer System according to the "dissection method" principle. Using singu- 
lar perturbation theory, we describe the sequence of bifurcations undergone by the regulating 
(FitzHugh-Nagumo) system, encompassing the rarely investigated case of homoclinic connex- 
ion. Basing on pure dynamical considerations, we restrict the space of parameter search for the 
regulating system and describe a foliation of this restricted space, whose leaves define constant 
duration ratios between the surge and the pulsatility phase in the whole system. We propose 
an algorithm to fix the parameter values to also meet the other prescribed ratios dealing with 
amplitude and frequency features of the secretion signal. We finally apply these results to 
illustrate the dynamics of GnRH secretion in the ovine species and the rhesus monkey. 
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1. Introduction 

In vertebrates, the reproductive system is made up of the hypothalamus, belonging to the 
central nervous system, the pituitary gland and the gonads (ovaries in females, testes in males). 
Within the hypothalamus, specific neurons secrete the GnRH (Gonadotrophin Releasing Hor- 
mone) in a pulsatile manner. This pulsatility has a fundamental role in the differential control of 
the secretion of both gonadotropins by the pituitary gland: LH (luteinizing hormone) and FSH 
(follicle stimulating hormone). In females, the pulsatile pattern is tremendously altered once per 
ovarian cycle into a surge which triggers LH surge and ovulation in response to increasing levels 
of estradiol secreted by the ovaries. The estradiol signal is conveyed to GnRH neurons through 
a network of interneurons. The balance between stimulatory and inhibitory signals emanating 
from interneurons controls the behavior of the GnRH network. 

We have proposed a mathematical model accounting for the alternating pulse and surge pattern 
of GnRH secretion [SJ. The model is based on the coupling between two systems running on 
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different time scaies. We thus consider the following four-dimensional dynamical system: 

(1.1a) eSx = -y + f(x) 

(1.1b) ey = a<jx + a\y + a 2 + cX 

(1.1c) ejX = -Y + g(X) 

(Lid) Y = b X + b 1 Y + b 2 



where / and g are two cubic functions : 

f(x) — v x 3 + uix 2 + v 2 x 
g(X) = fioX 3 + ^X 2 + fi 2 X 

The faster system (|l.la[) - (|l.lb[) . henceforth named "GnRH Secreting System", corresponds to 
the average activity of GnRH neurons, while the slower one (|Llc|l - (|l.ld|) . named "Regulating 
System", corresponds to the average activity of regulatory neurons. The x, X variables represent 
the neuron electrical activities (action potential), while the y,Y variables relate to ionic and 
secretory dynamics. In each system, the fast and slow variables feedback on each other. The 
coupling between both systems is mediated through the unilateral influence of the slow regulatory 
interneurons onto the fast GnRH ones (cX term in equation (|l.lb[) V The coupling term aggregates 
the global balance between inhibitory and stimulatory neuronal inputs onto the GnRH neurons. 
The analysis of the slow/fast dynamics exhibited within and between both systems allows to 
explain the different patterns (slow oscillations, fast oscillations and periodical surge) of GnRH 
secretion. 

This model can be used as a basis to understand the control exerted by ovarian steroids (estra- 
diol and progesterone) on GnRH secretion, in terms of amplitude and frequency of oscillations 
and discriminate a direct (on the GnRH network) from an indirect action (on the regulatory 
network) of steroid feedbacks. To account accurately for this control, we have to fully under- 
stand the sequence of bifurcations corresponding to the different phases of GnRH secretion and 
the occurence of an hysteresis loop. Our main goal in this paper is thus to carry on with a 
deeper analysis of this model, based on singular perturbation theory, to meet precise quantitative 
neuroendocrinological specifications. These specifications deal with the duration of the luteal 
(progesterone-dominated) and follicular (estradiol-dominated) phases of the ovarian cycle, and 
the ratios between (i) the surge duration and the whole cycle duration, (ii) the pulse amplitude 
and the surge amplitude, and (iii) the increase in pulse frequency from the luteal to the follicular 
phase. 

The paper is organized as follows. In section [2l we remind the main properties of the original 
model and explain the mechanisms underlying the alternation between pulsatility phases and 
surges of the secretion signal. In scction[31 we give a precise bifurcation diagram of the Regulating 
System according to the parameters e, b\ and b 2 , allowing to restrict the domain of parameter 
values. We prove the existence and describe the shape of a foliation of the restricted domain, 
whose leaves correspond to constant duration ratios. This novel type of control form enables us to 
find the parameter values fulfilling any such prescribed ratio. Section [5] deals with the influence of 
other parameters on the remaining specifications. It calls to methodological developments falling 
into the field of weakly coupled nonlinear oscillators. Finally, in section [SJ we enunce in details 
the quantitative specifications for the ovarian cycle in the sheep and rhesus monkey, since direct 
measurements of GnRH are available in these species. In each case, we instance the search for 
the relevant parameter values and illustrate the resulting GnRH secretion signal from numerical 
simulations. 
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2. Qualitative understanding of the original model 

In this section, we remind the main mechanisms underlying the model behavior, which is 
necessary to introduce the subsequent analysis in section 3. Following [3], we assume that all 
parameter values are positive except v$ and /xq- 

2.1. Reparameterization. Some parameters of (jl.ip are useless since, whatever their values, 
we can fix the other parameter values to obtain the same orbit. We can first remove the bo 
parameter through a simple rescaling, which transforms system (jl.ljl into: 

(2.1a) e'S'x = -y + f(x) 

(2.1b) e'y — a x + a^y + a 2 + cX 

(2.1c) e'jX = -Y + g{X) 

(2. Id) Y = X + b[Y + b' 2 

with: 

e = b s, b 1 = —, b 2 = — . 

We fix bo = 1 from now on. We can also remove the 7 parameter since system (jl.l[) also reads: 

(2.2a) e"5"x = -y + f(x) 

(2.2b) e"y = a'oX + a"y + a! 2 ' + c" X 

(2.2c) e"X = -Y + g(X) 

(2. 2d) Y = X + b 1 Y + b 2 



with: 



£7, 5" = -, a'o = a 7, a" = a-rf, a 2 = a 2 7, c" = cry. 



1 

We will fix 7 = 1 in the sequel. 

2.2. Qualitative behavior of the uncoupled system. The average activity of neurons of 
both regulating and secreting neurons is modeled by a classical FitzHugh-Nagumo system, while 
the control of GnRH neurons activity by regulatory neurons is expressed by means of a one-way 
coupling of variable X onto variable y. Hence, we can analyze separately the slow-fast subsystem 
(|l.lc[) - (|l.ld|) . according to the time scale parameter e and the F-nullcline parameters 60, bi, and 
b 2 . 

The classical approach of slow-fast systems consists in studying the "Boundary-Layer System" 
obtained after fixing the slow dynamics at (here Y = 0) and focusing on the geometric invariants 
of the fast dynamics. Precisely, with the time rescaling: 

(2.3) t = sr 

subsystem (fT7Ic]) - (|l.ld[) and (RS e ): 

(2.4a) X = -Y + g(X) 

(2.4b) Y = e(X + bxY + b 2 ) 

are conjugate. Posing e = in this last system, we obtain the corresponding Boundary-Layer 
System: 

(2.5a) X = -Y + g(X) 

(2.5b) Y = 
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This dynamics is a zero-order approximation of system (ll.lc|) - (|l.ldj) according to e. The singular 
perturbation theory allows us to assess where this approximation is valid. 

We consider the set of singular points of subsystem (|2.5p , defined by the cubic function Y = 
g{X), and analyze their nature according to the fast dynamics (|2.5a|) . The two local extrema 
of the cubic are non hyperbolic points. The points on the cubic middle branch (connecting 
the two local extrema) are hyperbolic repulsive, while the two other (left and right) branches 
consist of hyperbolic attractive points. Each branch is a normally invariant hyperbolic manifold 
of subsystem (12.51) . They persist under perturbation by the slow dynamics (i.e. for small values 
of e) into a 0(e)-close normally invariant hyperbolic manifold of (|l.lcp - (|l.ld[) . The attractive 
manifolds attract the current point, as soon as it lies within their attraction basin, so that the 
orbit actually enters a 0(exp(—k/e) neighborhood. 

In the classical case we are interested in, the values of parameters b\ and 62 are small enough 
for the y-nullcline to intersect the cubic on the middle branch at a repulsive singular point of 
system (| 1 . lc[) - (|1 . ld[) and the two other intersection points to be far away, respectively on the left 
and right branch. It is well-known that, in such a case, subsystem (|2.5[) admits an attractive 
limit cycle surrounding the middle singular point (cf Figure [TJa). 
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Figure 1. Instance of orbit projections and signals generated by (jTTTJ) . Each 
color corresponds to the same time mark in the four panels for which the same 
parameterization has been used, a) Attractive limit cycle of the Regulating Sys- 
tem (|l.la[) - (|l.lbp . b) X(t) output signal of the Regulating System, c) Projection 
of system orbit on the (x,y)-plane. d) y(t) output signal of the GnRH Se- 
creting system. 

The pulsatility phase duration (resp. surge duration) corresponds, within a 
0(e)-approximation, to the time during which X < (resp. X > 0) along the 
parameterization of the Regulating System limit cycle. 
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In the course of a period, the current point first comes near to, and then goes down slowly 
along the left perturbed manifold (near the left branch of the cubic Y — g(X)), at a 0(1) speed. 
Near the extremum of the cubic (which is a non hyperbolic singular point of (|2.5p ). it keeps 
on going down and moves away from the cubic. At that time, it is subject to the fast motion 
and quickly reaches (at a 0(1/ e) speed) the vicinity of the right branch, as its Incoordinate 
remains almost constant. Then it slowly goes up along the right perturbed manifold, before it 
finally comes back quickly near the left branch. This well-known motion is detailed in the proof 
of Lemma 13.11 For more details on singular perturbation theory, we refer to [5] for the general 
study of slow-fast systems and [T3] for its application to oscillations. 

2.3. Qualitative behavior of the whole system. We briefly remind here the specific behavior 
of system (jl.ip and the corresponding pattern of the y(t) signal generated as a model output. 
We refer to [3] for more explanations about the model properties. 

Through the coupling of (|1 . lajl - (|1 . lb|) with (|1 . lcjl - (| 1 . lcijl . the Regulating System drives the 
dynamics of the GnRH Secreting System, by controling the position of its straight-line separatrix 
defined by: 

aox + diy + a,2 + cX = 
When all parameter values are fixed, the dynamics is determined by X values. Let us consider 
X as a parameter. As its value increases, system (|l.la|) - (|l.lL)[) displays the following sequence of 
behaviors: 

(1) There is a unique singular point far up on the left branch. 

(2) An inverse saddle-node bifurcation occurs. 

(3) There are three singular points: a saddle far up on the left branch, an attractive point 
on the right branch and another saddle on its right. 

(4) A supercritical Hopf bifurcation occurs. 

(5) One of the three singular points lies on the middle branch and is surrounded by an 
attractive relaxation limit cycle. 

(6) An inverse supercritical Hopf bifurcation occurs. 

(7) There are three singular points: a saddle far down on the right branch, an attractive 
point on the left branch and another saddle on its left. Any orbit starting near the origin 
reaches the vicinity of the left branch and approaches the attractive singular point. 

(8) A saddle-node bifurcation occurs. 

(9) There is a unique singular point far down on the right branch. Any orbit starting near 
the origin reaches the vicinity of the left branch and goes up along. 

In the sequel, we will set precise conditions on the parameters such that, in the range of values 
taken by X along the limit cycle of the Regulating System, the GnRH Secreting System displays 
periodically the sub-sequence of behaviors described above from step 5 to 9 and back from step 
9 to 5. In this case, along an orbit of (jl.ip . the current point (x,y, X,Y) displays the following 
motion (represented in the panels of Figure [2] 

When (X,Y) goes down near the left branch of the cubic Y = g{X), X increases slowly: 
system (|l.laj) - (|l.lb|) displays the step-5 behavior and (x, y) turns around the limit cycle. Then, 
X increases very quickly, the limit cycle of system (|l.lal) - (jl.lb[) disappears (step 6). X decreases 
slowly and (x, y) climbs up the left branch of the cubic y = f(x) (step 7 to 9). Afterwards X 
decreases very quickly and (x, y) goes down near the left branch before oscillating again around 
the relaxation limit cycle. Typical orbits of system (II. ip generate y(t) signal patterns such as 
that represented in Figure 03-d. 

2.4. Pulse amplitude. From the previous discussion, we have seen that a pulsatile y(t) signal 
is generated from the oscillations of (x, y) around the limit cycle of the GnRH Secreting System. 
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Figure 2. GnRH Secreting System behavior as (X,Y) describes the limit cycle 
of the Regulating System. We have set Y — est (e zero-order approximation). 
The objects relating to the Regulating System are drawn in red or orange, those 
relating to the GnRH Secreting System are drawn in blue or green: the X- 
nullcline Y = g(X) is the orange cubic, the i-nullcline y — f(x) is the blue 
cubic, the F-nullcline is the red line, the y-nullcline is the green line, the red 
square represents the (X, Y) current point of the Regulating System, the green 
diamond represents the (x, y) current point of the GnRH Secreting System, 
a) Beginning of the pulsatility phase : the GnRH Secreting system displays 
the step 5) behavior of the sequence described in section 12.31 b) Middle of the 
pulsatility phase, c) Pulsatility to surge transition : the GnRH Secreting system 
displays the step 9) behavior, d) Surge rise, e) Surge decline : the GnRH 
Secreting system has come back to step 5). 

For all values of X leading to such oscillations, the shape of the limit cycle of system (|l.la[) - (|l.lbl) 
remains almost the same (see Figure [TJc)). Consequently, the pulse amplitude is approximatively 
equal to the difference between the y coordinates of the two extrema PI and p[_ of y = g(x). 
Thus, we can keep a unique relevant parameter to define the cubic function g and represent the 
amplitude. This reparameterization is described in more details in [3]. 
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To simplify the calculations and keep similar expressions for / and g, we pose: 

f[x) = -x 3 + 3X 2 x 
g(X) = -X 3 + 3n 2 X 

Then, the coordinates of the local extrema, P£ and P±, as well as those of their projections on 
the other branch of the cubics (along the trajectories of the fast systems), Q ± and Q±, are easy 
to compute: 

(2.6a) Pjf = (A, 2A 3 ), P[ = (-A, -2A 3 ) 

(2.6b) i* = (fx, 2fi 3 ), P 9 _ = (- M , -2//) 



(2.7a) Q f + = (-2A, 2A 3 ), Q f _ = (2A, -2A 3 ) 

(2.7b) Q\ - (-2/x, 2/. 3 ), Q 9 _ = (2/i, -2^ 3 ) 

These reference points and their coordinates are displayed in Figure [3] 




Figure 3. Graph of the cubics defining the x and X nullclincs. The coordinates 
of the local extrema P± (resp. P^) of the x-nullcline (resp. X-nullcline) are given 
in (|2.6[) . The projections Q£ (resp. Q±) of these extrema along y = est (resp. 
Y = est) are given in (|2.7p . 

Henceforth, we set A and fj, to given, constant values. We also fix S to a value small enough 
for the GnRH Secreting System to admit a relaxation limit cycle in the step-5 case described in 
section 12.31 Each of the other parameters constitute a degree of freedom, whose value has to be 
fitted to any specific list of specifications, as it is illustrated in ^51 We begin in the next section 
by studying the effect of e, &i and b 2 on the duration ratio between the surge and the whole cycle. 
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2.5. Control of the two-way transition between the pulsatility phase and the surge 
by the dynamics of the regulating system. From the sequence of behaviors described in 
W2.31 we can establish a direct link between the pattern of the GnRH secretion signal and the 
values of X along two distinct parts of the Regulating System limit cycle. 

As stated in [3] , the time during which the current point [X, Y) remains in a neighborhood 
of the right branch of Y = g(X) along the limit cycle is 0(1) in e. In contrast, the transition 
time of the Regulating System from a neighborhood of Pi to a neighborhood of Q 9 + is 0(e). The 
same e-order applies respectively to the slow motion near the left branch of the cubic Y = g{X) 
and the fast motion between P£ and Q 9 _. 

In a precise zone of the parameter space defined in $5j we can ensure that, for — 2fx < X < —fi, 
the GnRH secretion signal has a pulsatile pattern (step 5), while, for /i < X < 2/i, it switches 
to a surge pattern (step 7 to 9). Once {X,Y) has reached the vicinity of Q 9 + , it takes a 0(s)- 
time to (x, y) to reach a neighborhood of Q_ and oscillate again around the limit cycle of the 
GnRH Secreting System. Consequently, if we neglect the 0(e)-duration of fast motions, we can 
approximate the duration of the pulsatility phase by the time during which X < and the surge 
duration by the time during which X > (see panels c and d of Figure [1} . 

3. Bifurcation Diagram of the Regulating System 

In this section, we first describe the sequence of bifurcations undergone by the Regulating 
System (RS £ ), defined by (|2.4p and the corresponding behaviors exhibited by the whole system 
(whether such behaviors are meaningful or not from the biological viewpoint) in delimited zones 
of the parameter space. Then, we restrict the parameter space to a domain where the whole 
system meets the prescribed dynamics, alternating between pulsatility phases and surges. 

We are interested in determining conditions under which the Regulating System admits an 
attractive limit cycle, further denoted by C(bi,b2, s). We have first to assume that system (12.41) 
admits a singular point (Xq,Yq = g (Xq)) such that \Xq\ < [i, and that e is small enough (in 
a sense that we will explain later in Lemma 13. Yet, even under those necessary conditions, 
two different types of bifurcation may make the limit cycle disappear: a supercritical Hopf bi- 
furcation and a homoclinic bifurcation. To explain the occurrence of these bifurcations, we have 
to remind that, for any b\ > 0, system (|2.4p admits two other singular points of saddle type 
(X_,Y_ =g{X-)), X_ < —fi and (X+.Y+ = g(X+)), X + > /x. 

3.1. Hopf Bifurcation. The Hopf bifurcation happens when the middle singular point (Xq,Yo) 
coincides with the local minimum of the cubic Y = g(X), i.e. when Xo = — /i. Then : 

-fj, + &i5(-M) +h = -£« + h (y - 3^ 3 ) + b 2 = 

Hence, the Hopf bifurcation occurs along the surface: 

(3.1) Hp : b 2 = hpibx) = fi + 2b lf i 3 

in the parameter space (e, &i, 62). 

The Hopf bifurcation gives birth to a small limit cycle which surrounds the middle singular 
point (X ,Y ) as (61,62) crosses transversally H p . Since (RS e ) is a slow-fast system, this limit 
cycle becomes quickly a big relaxation limit cycle when (e, 61, 62) moves away from Ti p . 

The occurrence of small limit cycles in such relaxation systems is known as "Canard phenom- 
enon" . The "Canard cycles" are attractive limit cycles which follow, against all expectations, a 
neighborhood of the fast dynamics manifold of repulsive points for a while. Some examples are 
displayed in Figure 5] 
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Figure 4. Examples of Canard cycles of the Regulating System (RS e ) for three 
different yet very close (to 1CT 12 ) values of e. 

3.2. Homoclinic Bifurcation. The homoclinic bifurcation occurs when the left saddle point 
{X-,YJ) "crosses" the limit cycle, leading to the intersection of its stable manifold with its 
unstable manifold. We can intuitively understand this bifurcation by considering the common 
limit periodic set Tq approached by the family of periodic orbits for the Hausdorff distanccQ, 
when e — * 0. This set is the compact union of the cubic parts from Q 9 _ to and from Q 9 + 
to Pf on one hand, and the fast trajectories ]P+,Q+[ and ]P 9 ,Q 9 _[ on the other hand; it is 
represented in Figure O 

Hence the e, zero-order approximation of the limit cycle allows us to develop the zero-order 
approximation of the homoclinic connexion: 



It is equivalent to state that Q+ = (Xq + , g(Xq + )^ = (— 2/i, 2^ 3 ) lies on the separatrix X + b\Y + 



To delimit the domain in the parameter space where the e zero-order approximation is valid, 
we have to take into account the Canard phenomenon occurring near the Hopf bifurcation. In 
particular, the occurrence of small cycles such as those represented in Figure 2] shows that, in the 
region of the parameter space close to the Hopf bifurcation surface H. pi JiP c does not approximate 
the homoclinic bifurcation condition. The occurrence of homoclinic bifurcations in such a case is 
an awkward problem which we will not tackle in this article. Consequently, we make the suitable 

1 The Hausdorff distance between two compacts K and K' of a metric space {E, d) is the smallest r > such 
that Vx S K,d(x,K') < r and V» £ K',d(y, K) < r. 



(X_,g(X_)) = (X Q+ ,g(X Q+ )) 
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Figure 5. a) Limit cycle of the Regulating System (RS £ ) for the following values 
of parameters: b\ — 0.1, 6 2 = 1, e = 0.1. b) Limit Periodic Set for 6i = 0.1, 
62 = 1: limit of the family (C(£>i, 6 2 , e)) ee ]o ;£o ] as e tends to according to the 
Hausdorff distance. 

assumptions to remain in a zone of parameter values for which this question is not bothering. To 
do so, we need the following lemma: 

Lemma 3.1. For each a > 0, for all (61,62) verifying: 

(3.3) 6 2 < n + 2hn 3 - a 
and: 

(3.4) 62 < 2fi - 26i^ 3 

there exists Eq > such that, for all e 6]0, £o[, the limit cycle of {2.$ exists and contains some 
points of [/i,+oo[xR. Moreover, for (61,62) fixed, the limit cycle C{b\,b2,e) lies in a 0(e 2 ^ 3 )- 
neighborhood of Tq . 

Proof. We use the sketch of the well-known proof of limit cycle existence for the FitzHugh- 
Nagumo system, when the left and right singular points are respectively far away from Pf and 
Pf . This proof is based on the series expansion of the first return map induced by the flow of 
(RS S ). We focus here on the well-definition of this map. We base on the arguments enunciated 
in [T] to derive the series expansion of the map and prove the existence of one stable fixed point 
and that of the attractive limit cycle. The geometric objects mentioned below are represented in 
Figure El 

Let a > and (b%, 62) verifying (|3.3p and (|3.4|) . The part of the cubic between (X_, YL) and 
Pf, denoted by A4°_, is a normally attractive manifold of the boundary-layer system associated 
to (RS e ), which perturbs into a left 0(e)-close normally attractive manifold of (RS £ ), denoted 
by Ai s _. Let us consider a small section S transverse to M°_. Then, all trajectories of (RS e ) 
starting from E must remain in an exponentially small neighborhood of .ML until they enter a 
0(e 2/,3 )-neighborhood of Pf (see According to (|3.3p . the middle singular point of (RS E ) 

is above Q 9 ,. Thus, we can choose e small enough so that this neighborhood does not contain 
(Xo,Y ). 

Let us remind too that the limit cycle crosses X = —fj, at a point Pf below the point Pf , 
since y < along this part of the limit cycle. All trajectories of (RS e ) starting from a 0(e 2 / 3 )- 
neighborhood of Pf keep away from any singular point for e small enough. Hence, they remain in 
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Figure 6. The part M°_ of the cubic Y = g(X) between and P 9 _ is 

a normally invariant attractive manifold for the system (RSq). The part of the 
cubic Y = g(X) between Pf and (X + ,Y + ) (not shown on the picture, far down 
on the right branch) is a normally invariant attractive manifold for the system 
(RSq). Thus, A4°_ (resp. M+) perturbs into a 0(e)-close manifold M e _ (resp. 
.Ml), which is invariant and normally hyperbolic for (RS £ ). 
Starting from an initial data on section E transverse to Ai°_, we can track the 
orbit of (RS e ) until X — — //, then near M. E + until X = fx, and, at least, to E. 

a 0(e) -neighborhood of the fast trajectory starting from Pf , which reaches the right attractive 
manifold. 

With the same arguments as above, we keep on tracking the orbits of (RS e ) in an exponentially 
small neighborhood of the right 0(e)-close normally attractive manifold of (RS e ), denoted by 
A4 e + , until they reach a point Pf lying on X = // and above Pi. According to the theory 
of regular perturbations, the orbits which have just followed the perturbed attractive manifold 
remain in a 0(e 2 / 3 )-neighborhood of the fast trajectory connecting P^ with Q 3 + . From (|3.4[) . the 
left saddle point is strictly above the point Q+ on the cubic. This allows us to choose e small 
enough to make sure that all trajectories cross the cubic and thereafter remain in an exponentially 
small neighborhood of .ML until it crosses E. Thus the first return map is well-defined on E for 
e small enough, i.e. e €]0,£o[- The limit cycle C(bi,b2,e) lies in a 0(e 2 / 3 )-neighborhood of the 
limit periodic set Tq (see Figure [S]) and contains some points of [fj,, +00 [xM. □ 

Remark 3.1. The smaller a is, the smaller £q has to be, to avoid any Canard effect. Hence, 
from now on, a is fixed to a constant value small enough so that the region: 
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is non empty. 

On the other hand, the closer 62 is from 2/i — 26i/i 3 , the closer the left singular point is from 
the limit cycle, and the smaller Eq has to be. Eq tends to as (61,62) approaches Ti° c . 

In other words, for a fixed value of e, a homoclinic connexion of the left singular point occurs 
as (61,62) approaches 7i®. We seek this connexion as 62 increases from 0, assuming that: 

n -L n, 

6l > 



4/i 3 



When e is fixed to a small value, there exists (61,62), verifying (|3 . 3[) and (|3.4p . for which the 
value of £0, defined in lemma l3~Tl is greater than e. Then, for 62 small enough, the left singular 
point lies outside the 0(e 2 / 3 )-neighborhood of Tq containing the limit cycle C(6i,62,e). For a 
value of 62 close to (and smaller than) 2/i — 26i/i 3 , the saddle undergoes a homoclinic connexion 
where the limit cycle disappear. 

The following proposition enunciates the condition for a homoclinic bifurcation to occur away 
from any Canard situation : 

Proposition 3.1. There exists, locally near e = 0, a C 1 -surface in (61,62,6) of homoclinic 
connexions given by the graph: 



defined for: 



H c :b 2 = h c (bi,s) = 2/i - 2/x 3 6i 
jU + a 1 



O 



61 € 



4/i 3 V 



Proof. We have seen that the homoclinic connexion occurs when the left singular point, (X_, YJ), 
with X_ < —ji, crosses the limit cycle. Let X mnl (b±, 62, e) be the minimum of X along the limit 
cycle C(6i, 62, e). From the fast dynamics (which is null on the cubic), we have immediately that 
(^min(6i, 62, s), g(X m in{bi, 62, s)) lies on the cubic. Consequently, the condition for the homoclinic 
connexion reads: 

X min (h,b2,e) + b 1 g(X min (b 1 ,b 2 ,e)) + 6 2 = 
As we know that X min (6i, b 2 , e) = X Q+ + O (e 2 / 3 ) = -2/i + O (e 2 ^ 3 ), it also reads: 



(3.5) 



u(b x , 6 2 , e) = -2 M + 0(e 2/3 ) + 6i.g(-2/i + 0{e z/i )) + b 



2/3^ 







where u is a C 1 function on its definition domain (where the limit cycle exists) . For each 61 value 
such that: 



(3.6) 



61 € 



/i + a 1 



-Xmin(6i, 2/i — 26i/i 3 , 0) is well defined (it is equal to Xq + ) and: 

(ill 

— (6 1( 2 M -26i/i 3 ,0) = l 
<?6 2 

Then, the implicit function theorem implies the existence of a unique root of u with respect to 
62 for any 61 verifying (|3.6p and e small enough. From (|3.5p . this root can be expanded as: 



6 2 = 2/t - 2/i 3 6i + O (e 2/3 ) , e -> 



□ 
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3.3. Parameter space reduction. From the bifurcation diagram, we can reduce the parameter 
space in order to impose a periodic behaviour to the Regulating System. We can also ensure that, 
along its attractive periodic orbit, the variable X alternatively takes negative and positive values. 
Since the sign of X along the orbit of (|2.4p discriminates between the two phases (pulsatility phase 
and surge) of the whole system (jl.ip . the restriction of the value domain for (e, b\, b 2 ) is the first 
step to obtain an alternation of pulsatility phases and surges. 

In the sequel, we will only consider (£,61,62) to be under both the surface of homoclinic 
connexion and a security plane on the underside of the Hopf bifurcation: 



Tj — a 



b 2 = h p {bij 



Hence we assume that: 



{bi,b 2 ,e) e1li = { (b L . h 2 .:) 
where £0 is defined in Lemma T3. II 



b 2 < h p (b\) - a 
b 2 < h c (bi,e) 
e < e 




Figure 7. Bifurcation diagram for e < Eq- We restrict the parameter space to 
IZi : under the surface of homoclinic bifurcations TL C and far from the surface of 
Hopf bifurcation Tip. There, we are ensured that the Regulating System exhibits 
an attractive limit cycle and that, along this orbit, X takes alternatively positive 
and negative values. Both conditions are necessary to obtain the alternance of 
surges and pulsatility phases for the y(t) signal generated by the whole system 

tnj. 

In the next section, we will prove the possibility to tune the parameters 61, b 2 and e in this 
domain to fulfill the specification prescribing a given ratio of the duration of the pulsatility phase 
to the surge duration. 

4. Foliation of the Regulating System parameter space 

We focus in this section on the way of determining all the values of the triplet (61, 62, e) leading 
to a given ratio between the duration of the pulsatility phase and the surge duration. As we have 
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seen in §2.5[ these durations are related to the time during which the current point (X, Y) lies 
either in the half-plane X > (pulsatility phase) or X < (surge) within one period along the 
limit cycle of (RS £ ). Their ratio remains unchanged after time rescaling. 

Lemma 4.1. For (bi,b 2 ,e) £ TZ\, let us consider a parameterization (X(t),Y(t)) of the limit 
cycle C(&i,&2) s) such that X(0) = and Y(0) > 0, and let T(bi,b2,£) be its period. Then there 
exists a unique T_(&i,&2,e) such that: 

Vt e [0,T_(6 1 ,6 2 ,e)], X(t)<0 

Vt e [T_{b 1 ,b 2 ,e),T(b 1 ,b 2 ,e)}, X(t)>0 

Let us denote T + (6i, b 2 , e) = T(bi,b2,e) — T_(6i, b 2 , e). Then, as e tends to 0: 

t ( h k ^ r g'(X) + ^ r (X,e) , YMn( , 
r_(6i,o 2 ,e)= / T r ; — , — t 717 — ^ —dX + Oie) 

T + (bi,b 2 ,s)= Y , h( fyU , T , ,y rfX + Og 

where < cmd > are differ entiable functions defined respectively on: 

}X min (bi,b 2 ,s), -fi]x]0,e ) 



and: 



such that: 



[X min (bi,b2 ) e) ) -fi]x]0,e o ] 



(4.1) 3A(e) = o O(e 2 / 3 ),V1 e [X min (6i ) 6 2 ,e),-^] ) |* ± (X,e)| < A(e) 

Proof. We consider only the left part of the limit cycle, along which X < 0, and prove the 
existence of The argument for the right part is identical. 
Let us assume (b\, b 2 ,e) S IZi and consider the set: 

U = {{X, Y)\X <-fi,Y < g(X min (bi,b 2 , e), Y < g(X)} 

Since X = along the left branch of the critical manifold and Y < in U, any orbit of ()2.4I) 
starting in {J escapes from U across the half-line {(X,Y)\X = —fi,Y < 0}. From the dynamics 
(|2.4| . we have also, along any trajectory remaining in U : 



(4.2) ^ = y =£ x + &1 y + b 2 

V ' dX X -Y + g(X) 



It follows from 14.21 that such a trajectory can be represented as the graph of a differentiablc 
function of X. In particular, the limit cycle C(b 1 ,b 2 , e) goes through the point: 

(X m - m (bi , b 2 , e) , g(X min (6i , 6 2 , e) )) , 

enters U and remains in {/ until it crosses the half-line {(X, Y)\X = —fi, Y < 0}. Let us represent 
this trajectory as: 

C^(b 1 ,b 2 ,e) : Y = ^ but>2 ^(X), X e]X nlin (b 1 ,b 2 ,s), -fi] 

and pose: 

^:(X,s)^ 7bl , b2 , e (X)-g(X) 

We can directly deduce that \P_ < and is differentiable with respect to X. Since (bi,b 2 ,e) 
lies in IZi, the one-parameter family (C(bi, 62, s)) e e]o,e ] converges to the limit periodic set r 
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described in Lemma 13.11 and each C_(&i,&2,e) trajectory, e G]0, eo], remains in a O (e 2 / 3 )- 
neighborhood of {(X, g(X))\X €]X m i n (6i, b 2 , e), —£«]}■ Finally, , 5- is differentiable on: 

]X min (bi,b 2 ,e), -/i]x]0,e ] 

and verifies (|4.1[) . □ 

It is worth noticing that, as e tends to 0, the durations T_(6i,62,£) and X+ (61, 62, e) tend to the 
limit durations: 

(4-3) ^,62) = r—JtVL—dX 



(4.4) Tl{b u b 2 )= t 

J2, 



_ 2fl X + 6 l5 (X) + 6 2 



, a/1 X + 6 l5 (X) + 6 2 

which were chosen as approximations to compute the period of the limit cycle in [3], When 
(61, 62, e) is far under the surface of homoclinic bifurcations H c , it is indeed a good approximation 
(with a controlled 0(£ 2 / 3 )-error). But, in our case, since the duration of the pulsatility phase has 
to be much longer than the surge duration, we have to increase T_(6i, 6 2 , e) in comparison with 
T + (&i, &2> e)< To do so, the current point has to be confined for a while in the vicinity of the left 
singular point, where (bi, b 2 , e) may be very close to H. c . This is why we cannot neglect the time 
spent along the path from X = X m { n to X = — 2/x, even more so since the motion is very slow 
near the singular point. 

We now explain how to select the value of fe 2 from fixed values of 61 and e to meet any 
prescribed ratio = r. We first consider the case e = 0. 

Proposition 4.1. There exists a C 1 -foliation oflZi n {e = 0} of one dimensionnal leaves such 
that: 

1 ) each leaf is the graph: 

b 2 = l° r {h), 61 e k, \ 
L m 

0/ a differentiable function 1°, 

2) for each r > 1, </iere zs a Zea/ on which: 

(45) T^ M = 

Proof. From (14. 3|) and (]4.4[) . we can see that T° (61,62) and T?(6i,&2) are continuous functions 
in 7?.i n {e = 0} and differentiable with respect to (61, 62). Let us notice that for any value of 61: 

Tg(6i,0) 
r°(6i,0) 

This comes from the central symmetry of the dynamics with respect to the origin when 6 2 = 0. 

Besides, for any fixed value of 61, T° (61, 6 2 ) increases and T°(6i, 62) decreases as 62 increases. 
Thus: 

(4.6) «^W 1?(ll , W .«^M lS(lllW> o 

ab 2 db 2 

Let us first consider 61 fixed such that: 

a + u 1 

We can increase 62 until (61,62) reaches the line H c . When 62 tends to 2/j, — 2/z 3 6i, ^(61,62) 
tends to +00 whereas T°(6i,6 2 ) remains finite. Together with 14.61 this proves that, for each 
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r > 1, there exists a unique solution b 2 = lr(bl) m H {e = 0} of (|4.5|) . As T±(&i, 62) are 
C ll -functions with respect to (61, 62), the function Z° is C 1 with respect to i>x- 

On the other hand, let us consider the time ratio when (b±, 62) lies on the security plane H~ a , 
i.e.: 

a + \l 



bi < 



Then: 



b 2 = hpipx) - a 
9'(X) 



Tl(b 1 ,h p (b 1 )-a)= f 

J2 



-in x + big(X) + M + 2&1M 3 - a 
9'{X) 



2/i 



Let us define b\ as the unique solution in b\ of the equation: 

I?(6i,fcp(fti)-a) r 
Then, if we fix r > 1, for each b\ such as 

a + fi 



dX 
dX 



b\ < 61 < 



4^ 3 



there exists a unique solution 62 = of 14.5 



□ 



Theorem 4.1. There exists a C 1 -foliation oflZi, defined near e = 0, 0/ too dimensionnal leaves 
such that: 

1 ) each leaf is the graph: 



b 2 = l r {bi,e), 

of a differ entiahle function l r , 

2) for each r > 1, there is a leaf on which: 

T- (bi , b 2 , 



11 



T + (6i,62,e) 



Proof. We prolong T±(&i,&2)£) into a differentiable function at e = by posing T± (61, 62,0) 
7±(&i, 62). Let us fix r > 1 and b\ such that &^ < b\ < \. Then, like in the previous proof: 

ar_(5i,i?(6i),o)„ ru , 0/h , n , ar + (6 1) ij(6i),o)„ ;0 



-T + (6 1 ,/^(6 1 ),0) 



T_(6 1 ,^(6i),0)>0 



The implicit function theorem implies that there exists an open subset t/T i n the (e, &i)-space, 
e > 0, which contains (0,&i) and a differentiable function l r such that, for each (e, b\) in [7g : 

r_(6i,6a,e) 



admits a unique solution defined by: 



6 2 = f r (&l,e) 



By recovering all values of b\ such that b\ < b\ < -jy, we define Z r for 61 in 61, A 
enough. 



and eo small 

□ 
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For any prescribed ratio r between the duration of the pulsatility phase and the surge duration, 
we have thus shown that there exists a 2-dimcnsional manifold of solutions in the parameter space 
(61, 62, e), and we have provided a 0(e 2 / 3 ) approximation of the leaf by the plane: 




Figure 8. Leaves in (e, 61, 6 2 )-space defined by T_/T + = r, for r ~ 2,3, 6, 9, 12. 
This ratio corresponds to the pulsatility phase duration over the surge duration. 

We can find numerically the relevant surface solutions for different ratios, as it is illustrated 
in Figure [8] for r = 2,3,6,9 and 12. The highest r is, the closest the leaf is to the surface of 
homoclinic connexions. 

5. Tuning of the parameters of the GnRH Secreting System 

In this section, we describe the process for tuning the values of e, clq, ax, a%, and c with 
respect to a definite list of quantitative specifications. We first state the necessary conditions 
to ensure that system Ql.ip displays the alternation between pulsatility and surge described in 
i j2.3l Using geometrical tools, we then analyze the effect of each parameter on each prescribed 
ratio (duration of the surge over the whole cycle duration, pulse amplitude over surge amplitude, 
and increase in pulse frequency from the luteal to the follicular phase). Finally, we describe the 
algorithmic procedure to apply in order to tune the parameter values one after another and fulfill 
the specifications. 

5.1. Constraints on the parameters to keep the right order in the sequence of secre- 
tion patterns. We first study the parameter conditions under which the y(t) signal does start 
to pulse right after the surge. We have to prevent the GnRH Secreting System from exploring the 
step l-to-4 behaviors described in $2.3\ wherever (X, Y) is along the limit cycle of the Regulating- 
System. In other words, the Hopf bifurcation (step 4) of the GnRH Secretion System has to 
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occur for a value of X such that X < —2/i. This condition amounts to say that for X — — 2/i, y 

of 
+ ■ 



is positive at the extremum P(_, so that we have to impose the following inequality on c: 



a A + 2aiA 3 + a 2 

(5 ' 1} C " ^ 

We also have to ensure that the pulsatility phase does not break off before the Regulating 
System switches to the fast motion from P£ to Q 9 _ that induces the surge. The GnRH Secreting 
System has to keep on oscillating long enough, so that the step 6 reverse Hopf bifurcation cannot 
occur for values of X strictly smaller than — /x. This condition reads: 

-a A - 2aiA 3 + a 2 

(5.2) c > 

The GnRH Secreting System has to remain in the surge mode (where the y-nullcline lies on 
the right of P_ ) until X undergoes the fast motion between P£ and Q 9 + that brings the secreting 
system back to the pulsatile mode. This again imposes a condition on the reverse Hopf bifurcation 
that cannot occur for a value of X strictly greater than \x: 

, K x . a X + 2ai\ 3 - a 2 

(5.3) c > 

fi 

Bringing together the two last conditions, we obtain the following constraint: 

aoX + 2aiA 3 — a 2 



(5.4) c > 



Finally, we have to take care that the current point (x, y) does cross the y-nullcline from the 
left to the right when (X,Y) undergoes the fast motion from Pf to Q 9 + . To this end, the slope 
of the y-nullcline (a /ai) has to be steep enough, hence the value of a\ has to be small enough. 
Otherwise the current point (x, y) might escape from control at the end of the surge (step 9 
behavior described in q2.3[) and go on climbing up the left branch of the cubic y — f(x), on the 
left side of the y-nullcline. 



5.2. Effect of the e, ao and a 2 on the ratios prescribed to the secretion signal. 

Effect of e. When a\ is small enough, we can approximate the period of the GnRH Secreting 
System, for any value of X ranging between —2/i and —fi, by: 

(5.5, . ( r m — - dx + f m — ix \ 

\J_2\a x + aij{x)+a 2 + cX J 2X a x + aif(x) + a 2 + cX 



Under conditions (|5.ip and (|5.4p . this period is both well-defined and finite. We can deduce from 
expression (15. 5p that the pulse frequency increases as e decreases. Similarly, the surge amplitude 
increases as e decreases, since the (x, y) point moves along the left branch of the cubic y — f(x) 
at a 0(1/ e) speed, compared to the 0(l)-speed of the Y motion, that we use as reference. 
Effect of ao and a 2 . The surge amplitude increases exponentially as ao decreases. Thanks to 
this exponential dependency, it is easy to find an ao value compatible with condition (|5.ip and 
leading to a high-amplitude surge. a 2 controls the position of the stripe described by y in the 
(x, y)-plane, as X increases from —2/i to — /1. When a 2 increases, the stripe moves leftwards in 
the (x, y)-plane, so that both the pulse frequency and the surge amplitude increase. 
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Effect of c. c controls the width of this stripe. Since conditions (|5.ip and (|5.4j) reduce drastically 
the range of admissible values for c, they also lower its influence on the y dynamics along the left 
branch of the cubic y = f(x). Hence, changes in c do not impact much on the surge. In contrast, 
a small change in the width of the stripe impacts a lot on the set of limit cycles followed by (x, y) 
as X increases from — 2fi to — //. Finally, the greater c is (even subject to condition ()5 . 1|) ^ . the 
greater are the stripe width and the pulse frequency ratio between the beginning and the end of 
the pulsatility phase. 



5.3. Procedure for tuning the parameter values. As it is explained in H2.41 we have fixed 
once and for all the values of S, A and p. More precisely, 6 — 0.0125 and the values of A and p 
correspond to cubic functions / and g expressed as: 

(5.6a) f(x) = -x 3 + 2.5.x 

(5.6b) g{X) = -X 3 + AX 

We propose here an algorithm-like procedure for tuning the other parameters in order to meet 
the set of quantitative specifications together. The procedure consists of the following steps: 

(1) Fix the value of a 2l since this parameter has multiple influences on the different ratios. 
To ensure that conditions (|5.1[) and (|5.4p will be fulfilled, we have to position the y- 
nullclinc precisely in the (x, y) plan in the uncoupled case (X = 0), so that it separates 
the P_ extremum, on its left, from the intersection point between the x-axis and the 
cubic y = f(x), on its right. The corresponding a 2 values range between A and ^/(3)A. 
Within this interval, 02 will be even greater that the prescribed surge amplitude is high. 

(2) Choose the order of magnitude of e to fit the average pulse frequency. Assuming that 01 
is small, we may approximate the minimum period of the GnRH Secreting System, for 
X ranging between —2p and —fi, by: 

(5.7) T min = 2e [ ^-dx = — (9A 2 - 6 ln(2)) 



2A a x a 

Hence, to obtain a prescribed frequency <f) at the end of the pulsatility phase, we can link 
ao and e by: 

(5.8) T mm = 1/0 

(3) Tunc the value of ao, subject to condition (|5.8jl . to fit the surge amplitude to obtain the 
suitable surge amplitude. 

(4) Find the value of c consistent with the pulse frequency ratio between the beginning and 
the end of the pulsatile phase. With 02 ranging between A and a/(3)A, the period of 
the limit cycle of the GnRH Secreting System for X = — p can be approximated by the 
minimum T m i n , whatever the value of c subject to conditions (|5.ip and (|5.4|) . Then, the 
period of the limit cycle for X = —2fx is equal to T m i n / p, where p denotes the frequency 
ratio. Finding c thus amounts to solve the implicit equation: 

- d x + r x m *, - T - 



(5.9) 



_ 2A a x + aif(x) + a 2 ~ cp J 2X a x + a\f(x) + a 2 - cp p 

It is worth noticing that this equation does not admit a solution in c for every value of 
p. The greatest ratio can be reached with the following value of c: 

aoA + 2aiA 3 — a 2 
c = 
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(5) Define the precise value of a±. We have already assumed that a\ is small enough (in a 
sense detailed in §5.1|) . The precise choice of ai affects marginally the amplitude of the 
surge, which increases as a\ increases 

(6) Deduce the values of b\ and 62 from the results established in $4] For a prescribed 
duration ratio, there is a 1-dimensionnal curve of solutions in the (&i, &2)-space. Along 
one such curve of constant ratio, the smaller b\ is, the longer X remains close to — /z, in 
comparison with the time spent near X — — 2/i, hence the sooner the pulse frequency 
increases. This property is used to set approximately the durations of the luteal and 
follicular phases, even if it is not straightforward to determine the value of X (between 
—2/i and —pi) for which the pulse frequency starts increasing drastically. 

6. Quantitative neuroendocrinological specifications 
6.1. GnRH secretion along the ovarian cycle in the ovine species. 

6.1.1. Duration of ovarian cycle phases. During the breeding season (late summer to the start 
of the next spring), the ewe shows ovarian (estrous) cycles of 16-17 days [8|. The cycle is divided 
into 2 phases, the follicular and the luteal phase. The follicular phase has a duration of 2-3 
days and is characterized by increasing secretion of estradiol and onset of the LH (Luteinizing 
Hormone) surge (triggered by the GnRH surge). The luteal phase has a duration of 14-15 days 
and is characterized by the secretion of progesterone from the corpus luteum. The transition 
from the follicular to the luteal phase is marked by ovulation. In the absence of pregnancy, the 
transition from the luteal to the follicular phase is marked by the luteolysis of the corpus luteum, 
allowing resumption of a new cycle. 

6.1.2. Amplitude and frequency ratios. We first detail the specifications concerning the ovine 
species, since this is the species for which there are the more data available directly on the 
GnRH level (rather than the LH level). Indeed, in the ewe, the development of a dedicated 
technique has allowed the sampling of pituitary portal blood with high temporal resolution [10] . 
This technique is specially useful in studying the pattern of GnRH secretion during the surge. 
It has been utilized both during the luteal and follicular phase of naturally cycling animals 
[llj and during artificial, estradiol- induced, follicular phases in ovariectomized ewes |11[ 19], In 
any case, the LH surge induced by estradiol was invariably accompanied by a massive GnRH 
surge. Depending on studies and between-animal variability, the peak values of the GnRH surgj^ 
ranged from 85 to 170 (73 to 394 outside the breeding season) [TO], 9 to 109 [TTJ fold over 
baseline, with 10 min sampling intervals, or even from 100 to 500 fold over baseline with 2 min 
sampling intervals [9]. This variability can be partly explained by differences in the anatomical 
level where the surgical cut penetrates the portal vasculature (the higher the cut the higher the 
surge amplitude). Comparatively, the average pulse peak values can be assessed as less than 10 
fold the baseline. Accordingly, we chose a surge to pulse amplitude ratio instance of 60 in the 
numerical simulation illustrated below. These studies had also established that the GnRH pulse 
frequency not only was greater in the early follicular phase than in the luteal phase, but also 
increased further within the follicular phase as the surge approached. More precisely, during the 
luteal phase, the average frequency was one pulse per 4 hours, while it exceeded one pulse per 
hour in the follicular phase. Accordingly, we considered an average pulse frequency ratio of 4 
between the luteal and follicular phase in the model. 



"In the cited papers, the surge peak value is actually an average computed from the values observed during one 
hour after the maximum has occured. Hence, with a 10 min sampling interval, this is an average over 6 values. 
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6.1.3. Parameter combination for the ovine cycle. Gathering the information detailed above, we 
set ourselves the target of meeting the specifications detailed in Table [TJ Applying the procedure 
described in we deduce the set of parameters listed in Table [5J The corresponding GnRH 
secretion pattern is illustrated in Figure [9l The zoom around the surge in Figure [10] emphasizes 
the difference in the pulse frequency at the end of the follicular phase (pre-surge period) compared 
to the beginning of the luteal phase (post-durge period). 

Table 1 . Quantitative specifications for GnRH secretion pattern along the ovar- 
ian cycle in the ewe 



feature 


specification 


unit 


whole cycle duration 


16.5 


days 


follicular phase duration 


2.5 


days 


surge duration 


1 


days 


luteal phase duration 


13 


days 


pulse to surge amplitude ratio 


60 




frequency increase ratio 


4 





Table 2. Parameter values for the ovarian cycle in the ewe 



£ 


0.02 


S 


0.0125 




0.52 


oi 


0.011 


a-2 


1.14 


c 


0.70 


h 


0.246 


b 2 


1.5103 



120 ■ 



103 ■ 



y(t) 




t (scale : 1 day) 

Figure 9. Ovarian cycle in the ewe. y(t) signal generated by (|l.lj) with the 
parameter values given in Table [D The time scale unit is one day. The y 
variable has no dimension as we are interested in the pulse to surge amplitude 
ratio. The signal respects the specifications for the sheep ovarian cycle : whole 
cycle duration of 16.5 days, follicular phase of 2.5 days, surge duration of 1 day. 
The pulse to surge amplitude ratio is around 60. 
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y(t) * ■ 




t (scale : 1 day) 

Figure 10. Zoom around the surge (ewe). y(t) signal generated by (jl.ip with 
the parameter values given in Tabled] The whole follicular phase (2.5 days) is 
represented here. The pulse frequency before the surge is approximatively 1 per 
70 minutes. The pulse frequency at the beginning of the luteal phase is around 
1 per 2.5 hours. 



6.2. GnRH secretion along the ovarian cycle in the rhesus monkey. 

6.2.1. Specifications on cycle phase duration and GnRH surge to pulse ratios. The duration of 
the ovarian (menstrual) cycle in the rhesus monkey is comparable to that of the human cycle, 
i.e around 28 days long. Comparatively to the ovine cycle, it is much more symmetric, since the 
average duration of both the follicular and the luteal phase amounts to 14 days 

GnRH measurements in cerebrospinal fluid (CSF) samples have been obtained from the third 
ventricle of intact and ovariectomized conscious rhesus monkeys during control periods and 
throughout an estrogen challenge [14], with a time resolution of 15 min. In the intact as well as 
ovariectomized rhesus monkey, a genuine GnRH surge does occur in response to estradiol and 
the profile of the GnRH surge is remarkably similar to that reported in the ewe. These semi- 
quantitative information are quite reliable since simultaneous measurements of GnRH in CSF 
and portal blood in the ewe have shown that there is a good agreement between both techniques 
at the time of the GnRH surge [12] . However, it is more difficult to extract accurate quantitative 
information from these data, since the technique of collecting GnRH in the CSF is less reliable 
than the sampling of pituitary portal blood, due to difficulties in maintaining the required CSF 
flow uninterruptedly for a long time and documenting the position of the tip of the collecting 
cannula. Moreover, the whole surge duration was not entirely covered by the collecting period 
in several monkeys, so that maximal values may have been missed. Considering the maximal 
observed GnRH CSF concentration (around 400 pg/ml) and a median pulse amplitude of 16 
pg/ml, we can however fix the pulse to surge amplitude ratio for instance to 25. 

Table 3. Quantitative specifications for GnRH secretion pattern along the ovar- 
ian cycle in the rhesus monkey 



feature 


specification 


unit 


whole cycle duration 


28 


days 


follicular phase duration 


13 


days 


surge duration 


1 


days 


luteal phase duration 


14 


days 


pulse to surge amplitude ratio 


25 




frequency increase ratio 


4 
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6.2.2. Parameter combination for the rhesus monkey cycle. Gathering the information detailed 
above, we set ourselves the target of meeting the specifications detailed in Table [3] Applying the 
procedure described in 5J5] we deduce the set of parameters listed in Table SJ The corresponding 
GnRH secretion pattern is illustrated in Figure [11] One can notice the difference in the length 
of the follicular phase compared to that of the ovarian cycle in the ewe. The zoom around the 
surge in Figure [12] emphasizes the difference in the pulse frequency at the end of the follicular 
phase (pre-surge period) compared to the beginning of the luteal phase (post-durge period). 

Table 4. Parameter values for the ovarian cycle in the rhesus monkey 



£ 


0.018 


6 


0.0125 


ao 


0.7 


a\ 


0.013 




1.0 


c 


0.67 


h 


0.187 


b 2 


1.704 



y(t) 




t (scale : 1 day) 



Figure 11. Ovarian cycle in the rhesus monkey. y(t) signal generated by (jl.ip 
with the parameter values given in Table [4] The time scale unit is one day. The 
y variable has no dimension as we are interested in the pulse to surge amplitude 
ratio. The signal respects the specifications for the rhesus monkey ovarian cycle 
: whole cycle duration of 28 days, follicular phase of 13 days, surge duration of 
1 day. The pulse to surge amplitude ratio is around 25. 



y(t)> 















III 


11 





t (scale : 1 day} 

Figure 12. Zoom around the surge (rhesus monkey). y(t) signal generated by 
(jl.ip with the parameter values given in Table [4] The three last days of the 
follicular phase are represented here. The pulse frequency before the surge is 
approximatively 1 per 80 minutes. The pulse frequency at the beginning of the 
luteal phase is around 1 per 2 hours. 



24 



FREDERIQUE CLEMENT AND ALEXANDRE VIDAL 



7. Conclusion and discussion 

We have based our study on a concise model reproducing the different GnRH secretory patterns 
along an ovarian cycle [3J. In this model, the dynamical pattern of GnRH secretion results 
from the interaction between two neuronal populations: GnRH secreting neurons and regulating 
neurons, whose paragon can be embodied by Kiss-peptin neurons [4]. The latter population 
integrates much of the ovarian steroid feedback and acts as a slow pacemaker, while the former 
one is part time excitable, part time fast oscillating. The alternation of pulsatile and surge mode 
for GnRH secretion, as well as pulse frequency increase as the surge onset gets closer, have been 
explained in the framework of bifurcation theory. 

In this paper, we have conducted a deeper analysis of this model, with the definite aim of 
constraining the model outputs (and more precisely the output variable corresponding to GnRH 
secretion), with respect to a physiologically relevant list of specifications. In other words, we 
have challenged the ability of the model to meet precise quantitative relations on the secretion 
signal features. Apart from the total duration of the ovarian cycle, which is expressed in physical 
time, these relations can all be expressed as ratios, regarding (i) surge duration over the whole 
cycle duration, (ii) the duration of the luteal phase over that of the follicular phase, (iii) pulse 
amplitude over surge amplitude, and (iv) pulse frequency in luteal phase compared to follicular 
phase. 

We have described in great details the sequence of bifurcations undergone by the uncoupled 
FitzHugh-Nagumo subsystems, beyond the usually-investigated situation where the slope of the 
slow variable nullcline is steep. Using singular perturbation theory and dynamical analysis, we 
have formulated a e-expansion of the homoclinic bifurcation surface in the (&i, 62, e) space. From 
this expansion, we have been able to restrict the space of parameter values search precisely, 
without reducing the set of reachable ratios. Within this restricted space, we have described a 
foliation, whose leaves define constant duration ratios between the surge and whole cycle. From 
then on, we have used this foliation to determine a sufficient condition to fulfill the specification 
on such a ratio. This condition links together the values of 3 parameters (61, 62, e) over the 7 to 
be fixed in the adimensionned form of the model, that we have rewritten in the beginning of the 
article. The remaining parameters can be further tuned according to an algorithm taking into 
account the other prescribed ratios. We have managed to cope with all the ratios, even if the 
frequency increase is structurally limited in the framework of the model. 

We have finally applied our results to reproduce the GnRH secretion pattern in two different 
species in which GnRH data are both directly available and reliable. The ovine species exemplifies 
an estrous cycle, whose most obvious sign consists of heat behavior, while the rhesus monkey 
exemplifies a menstrual cycle, whose most obvious sign consists of menstruation. The former 
species is mostly interesting from an agronomic viewpoint, while the latter is interesting from a 
comparative physiology viewpoint, since its ovarian cycle is the closest to the human cycle. Even 
in species for which fewer GnRH data are available, we can make use of our parameter search 
algorithm. By instance, if we assume that the GnRH secretion dynamics in cows is comparable 
to that observed in ewes, as suggested in [B] , and adapt the durations of the follicular and luteal 
phases to values appropriate for the cow [7] , we can bend the ovine GnRH signal into a bovine- like 
signal. In the same spirit, we plan to represent various physiological and pathological situations 
from our deep understanding of the model behavior. 

From a dynamical viewpoint, this study indirectly addresses the questions of tracking homo- 
clinic connexions, where classical Canard cycles disappear, and Bogdanov-Takens like bifurcations 
of Relaxation Oscillators that occur on a parameter control manifold. Future work will also fo- 
cus on determining whether the whole system admits a strictly periodic attractive orbit. It is a 
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challenging question dealing with the synchronization of weakly coupled oscillators and delay to 
bifurcation analysis. 
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